#!/usr/bin/python3
##################################################
'''求解一元n次方程组'''
# Copyright (c) 2022 Xu Ruijun | 1687701765@qq.com
# the file is MIT License
##################################################
from math import *

__author__ = "Xu Ruijun"
__copyright__ = "Copyright (c) 2022 Xu Ruijun"
__license__ = "MIT"


def solve_x2(a, b, c):
    delta = b**2 - 4*a*c
    if delta < 0:
        raise ValueError('no real solve')
    sqd = sqrt(delta)
    return (-b + sqd)/(2*a), (-b - sqd)/(2*a)

def cbrt(x):
    if x < 0:
        return -cbrt(-x)
    return pow(x, 1/3)

def solve_x3_shenjing(a, b, c, d):
    A = b**2 - 3*a*c
    B = b*c - 9*a*d
    C = c**2 - 3*b*d
    delta = B**2 - 4*A*C
    if delta == 0:
        assert -b/(3*a) == -c/b == -3*d/c
        return -c/b
    elif delta > 0:
        Y1 = A*b + 3*a*(-B+sqrt(delta))/2
        Y2 = A*b + 3*a*(-B-sqrt(delta))/2
        return (-b-(cbrt(Y1)+cbrt(Y2)))/(3*a)
    elif delta < 0:
        T = (2*A*b - 3*a*B)/(2*sqrt(A**3))
        theta = acos(T)
        X1 = -b-2*sqrt(A)*cos(theta/3)/(3*a)
        X2 = -b+sqrt(A)*(cos(theta/3)+sqrt(3)*sin(theta/3))/(3*a)
        X3 = -b+sqrt(A)*(cos(theta/3)-sqrt(3)*sin(theta/3))/(3*a)
        return X1, X2, X3
    else:
        raise ValueError(f'delta={delta}')

def solve_x3_cardano(a, b, c, d):
    w = (-1+sqrt(3)*1j)/2
    w2 = w**2
    p = (3*a*c - b**2)/(3*a**2)
    q = (27*a**2*d - 9*a*b*c + 2*b**3)/(27*a**3)
    q2 = q/2
    p3 = p/3
    sqd = sqrt(q2**2 + p3**3)
    u = -cbrt(q2+sqd)
    v = -cbrt(q2-sqd)
    y1 =    u +    v
    y2 = w *u + w2*v
    y3 = w2*u + w *v
    x1 = y1 - b/(3*a)
    x2 = y2 - b/(3*a)
    x3 = y3 - b/(3*a)
    return x1, x2, x3
